import numpy as np


# This is problem 1
def gaussy(A, b, n):
    s = [0] * n
    for j in range(0, n):
        s[j] = abs(max(A[j], key=abs))
    l = [i for i in range(n)]
    for k in range(0, n):
        ratiomax = 0
        maxrow = k
        for i in range(k, n):
            ratio = abs(A[l[i]][k] / s[l[i]])
            if (ratio > ratiomax):
                ratiomax = ratio
                maxrow = i
        temp = l[maxrow]
        l[maxrow] = l[k]
        l[k] = temp
        for j in range(k + 1, n):
            multiplier = A[l[j]][k] / A[l[k]][k]
            for i in range(k + 1, n):
                A[l[j]][i] = A[l[j]][i] - (multiplier * A[l[k]][i])
            b[l[j]] = b[l[j]] - (multiplier * b[l[k]])
    #print('This is A', A)
    #print('This is b', b)
    x = [0] * n
    x[n - 1] = b[l[n - 1]] / A[l[n - 1]][n - 1]
    for j in range(2, n + 1):
        summ = 0.0
        for k in range(n - 1, n - j, -1):
            summ = summ + A[l[n - j]][k] * x[k]
        x[n - j] = (b[l[n - j]] - summ) / A[l[n - j]][n - j]
    print('The solution vector is', x)


gaussy([[2, 3, -4, 1], [1, -1, 0, -2], [3, 3, 4, 3], [4, 1, 0, 4]],
       [-15, 6, 6, -5], 4)
print('\n')

# This is problem 1a
gaussy([[3, -13, 9, 3], [-6, 4, 1, -18], [6, -2, 2, 4], [12, -8, 6, 10]],
       [-19, -34, 16, 26], 4)
print('\n')

# This is problem 1b
gaussy([[2, 4, -2], [1, 3, 4], [5, 2, 0]], [6, -1, 2], 3)
print('\n')

# This is problem 1c
# n = 6
n = 6
Aa = np.zeros((n, n))
for j in range(0, n):
    for i in range(0, n):
        Aa[i][j] = 1 / (i + j + 1)

bb = np.zeros((n, 1))
for i in range(0, n):
    bb[i] = sum(Aa[i])
gaussy(Aa, bb, 6)

# n = 15
n = 15
Aa = np.zeros((n, n))
for j in range(0, n):
    for i in range(0, n):
        Aa[i][j] = 1 / (i + j + 1)

bb = np.zeros((n, 1))
for i in range(0, n):
    bb[i] = sum(Aa[i])
gaussy(Aa, bb, 15)

# This makes sense because when n=15 the condition number is
# much larger than when n=6. This means there is much more room
# for potential error so it makes sense that the solution vector
# doesn't add up to all 1's.


# This is problem 2
def seidel(A, b, n, x0, tol):
    maxiteration = 200
    iter = 0.0
    abserror = tol + 1
    while (abserror > tol and iter < maxiteration):
        iter = iter + 1
        abserror = 0.0
        xold = list(x0)
        for i in range(0, n):
            summ = 0.0
            for j in range(0, i):
                summ += (A[i][j] * x0[j])
            for j in range(i + 1, n):
                summ += (A[i][j] * x0[j])
            x0[i] = (b[i] - summ) / A[i][i]
            abserror += abs(x0[i] - xold[i])
    if (abserror > tol):
        print('Did not converge')
    else:
        print('Converges to', x0)


seidel([[3, 1], [2, 6]], [5, 9], 2, [1, 1], 10**(-2))
seidel([[9, -3], [-2, 8]], [6, -4], 2, [0, 0], 10**(-2))
seidel([[1, -3], [-2, 8]], [6, -4], 2, [0, 0], 10**(-2))
# There is a relatively large change in the solutions when you
# change just one number. This could be because there's a lot of
# room for potential error. This could be due to a possible
# large condition number.